其他
非人物种的GSEA&GSVA分析
生信会客厅交流群里有几个朋友问我：MSigDBS数据库只有人的基因集gmt文件可供下载，小鼠的GSEA&GSVA分析如何做？解决问题的根本是做同源基因转换，代码能力强的朋友很容易实现，在线转换的网站也是有的。今天我给大家介绍一个R包——msigdbr，可以更轻松地帮大家得到部分非人物种的gmt文件或基因集列表。
演示数据
演示数据来自GSE109774中的一个小鼠心脏子数据集，整理好的seurat对象文件scRNAmm.RDS，我已上传到百度云。
链接：https://pan.baidu.com/s/1kWrL56O1EkAHPAY3RC4E2g
提取码：29pj
GSEA分析
准备gsea输入文件
##R语言准备gsea输入文件library(Seurat)library(tidyverse)dir.create("GSEA") dir.create("GSEA/input")dir.create("GSEA/output")scRNA <- readRDS("scRNAmm.RDS")#提取endothelial cell和fibroblast细胞做富集分析tmp <- scRNA@meta.datatmp <- subset(tmp, subset = (tmp$cell_ontology_class=='endothelial cell'|tmp$cell_ontology_class=='fibroblast'))sub.cells <- rownames(tmp)scRNAsub <- subset(scRNA, cells=sub.cells)scRNAsub$cell_ontology_class <- gsub('endothelial cell','endothelial_cell',scRNAsub$cell_ontology_class)expr <- GetAssayData(scRNAsub, slot = 'data')expr <- data.frame(NAME=rownames(expr), Description=rep('na', nrow(expr)), expr, stringsAsFactors=F)write('#1.2', "GSEA/input/expr.gct", ncolumns=1)write(c(nrow(expr),(ncol(expr)-2)), "GSEA/input/expr.gct", ncolumns=2, append=T, sep='\t')write.table(expr, "GSEA/input/expr.gct", row.names=F, sep='\t', append=T, quote=F)line.1 <- c((ncol(expr)-2), 2, 1)tmp <- table(as.character(scRNAsub@meta.data$cell_ontology_class))line.2 <- c("#", names(tmp))line.3 <- c(rep(names(tmp)[1],tmp[1]), rep(names(tmp)[2],tmp[2]))write(line.1, 'GSEA/input/group.cls', ncolumns=length(line.1), append=T, sep='\t')write(line.2, 'GSEA/input/group.cls', ncolumns=length(line.2), append=T, sep='\t')write(line.3, 'GSEA/input/group.cls', ncolumns=length(line.3), append=T, sep='\t')准备小鼠gmt文件
##安装msigdbr包install.packages("msigdbr")#查看msigdbr包支持的物种msigdbr_species() species_name species_common_name 1 Bos taurus cattle 2 Caenorhabditis elegans roundworm 3 Canis lupus familiaris dog 4 Danio rerio zebrafish 5 Drosophila melanogaster fruit fly 6 Gallus gallus chicken 7 Homo sapiens human 8 Mus musculus house mouse 9 Rattus norvegicus Norway rat 10 Saccharomyces cerevisiae baker's or brewer's yeast11 Sus scrofa pig
##准备小鼠gmt文件，以GO数据库的CC数据集为例library(msigdbr)#选择自己需要的基因集，category对应MSigDB的主分类，subcategory对应次级分类genesets = msigdbr(species = "Mus musculus", category = "C5", subcategory = "CC") #整理成gmt文件do_gmt <- select(genesets, c("gs_name","gs_url","gene_symbol")) %>% as.data.frame()gs_url <- do_gmt[,c("gs_name","gs_url")] %>% unique.data.frame()do_gmt <- split(do_gmt$gene_symbol, do_gmt$gs_name)for(i in seq_along(do_gmt)){ c1 <- names(do_gmt)[i] c2 <- gs_url[which(gs_url$gs_name==c1),2] c3 <- do_gmt[[i]] cb <- c(c1,c2,c3) write(cb, 'GSEA/input/genesets.gmt', ncolumns=length(cb), append=T, sep='\t')}推荐使用桌面版GSEA软件分析，使用方法参阅《单细胞转录组高级分析五：GSEA与GSVA分析》，部分结果如下：
GSVA分析
使用msigdbr包创建gsea分析使用的gmt文件还有点麻烦，但是用它创建gsva分析使用的基因集列表非常方便。
library(Seurat)library(GSVA)library(msigdbr)library(tidyverse)dir.create("GSVA")##选择自己需要的基因集genesets <- msigdbr(species = "Mus musculus", category = "C8") %>% select("gs_name","gene_symbol") %>% as.data.frame()genesets <- split(genesets$gene_symbol, genesets$gs_name)
##提取表达矩阵scRNA <- readRDS("scRNAmm.RDS")expr <- as.matrix(scRNA@assays$RNA@data)expr <- expr[rowSums(expr)>=1,] #选取非零基因
##gsva默认开启全部线程计算，建议count用"Poisson"，data用"Gaussian"es.matrix = gsva(expr, genesets, method="ssgsea", kcdf="Gaussian", parallel.sz=20) write.table(es.matrix, 'GSVA/gsva.xls', row.names=T, col.names=NA, sep='\t')
##差异分析#挑选差异比较明显的genesets，过程略write.table(s.gsva, 'GSVA/s_gsva.xls', row.names=T, col.names=NA, sep='\t')
##热图可视化library(pheatmap)s.gsva <- read.table('GSVA/s_gsva.xls', header = T, row.names = 1, sep = '\t', check.names = F)pheatmap(s.gsva, show_rownames=1, show_colnames=0, fontsize_row=8, filename='GSVA/gsva_demo.png', width=15, height=10)我特意挑选了MSigDB数据库的C8数据集——C8: cell type signature gene sets：Gene sets that contain curated cluster markers for cell types identified in single-cell sequencing studies of human tissue. 感觉它能辅助鉴定细胞类型，从我的结果来看不太靠谱。
交流探讨
如果您阅读此文有所疑惑，或有不同见解，亦或其他问题，欢迎添加我的微信探讨。我工作的公司实力雄厚，拥有10年测序服务经验，欢迎大家联系Kinesin洽谈测序及数据分析业务！